library(tidyverse)
library(triangle)
library(patchwork)
library(plotly)
plot_shaded_normal <- function(mean, sd, p = 0.5) {
shade_up_to_x <- qnorm(p, mean = mean, sd = sd)
data_normal <- data.frame(
x = seq(mean - 4 * sd, mean + 4 * sd, length.out = 300)
) |>
mutate(y = dnorm(x, mean = mean, sd = sd))
shaded_data_normal <- data_normal |>
filter(x <= shade_up_to_x)
ggplot(data_normal, aes(x = x, y = y)) +
geom_line() +
geom_area(
data = shaded_data_normal,
aes(x = x, y = y),
fill = "#ff8800",
alpha = 0.3
) +
geom_vline(xintercept = shade_up_to_x)
}
plot_shaded_triangle <- function(a, b, c = (a + b) / 2, p = 0.5) {
shade_up_to_x <- triangle::qtriangle(p, a = a, b = b, c = c)
data_tri <- data.frame(
x = seq(a, b, length.out = 300)
) |>
mutate(y = triangle::dtriangle(x, a = a, b = b, c = c))
shaded_data_tri <- data_tri |>
filter(x <= shade_up_to_x)
ggplot(data_tri, aes(x = x, y = y)) +
geom_line() +
geom_area(
data = shaded_data_tri,
aes(x = x, y = y),
fill = "#ff8800",
alpha = 0.3
) +
geom_vline(xintercept = shade_up_to_x)
}
plot_shaded_unif <- function(a = 0, b = 1, p = 0.5) {
shade_up_to_x <- qunif(p, min = a, max = b)
data_unif <- data.frame(
x = seq(a, b, length.out = 300)
) |>
mutate(y = dunif(x, min = a, max = b))
shaded_data_unif <- data_unif |>
filter(x <= shade_up_to_x)
ggplot(data_unif, aes(x = x, y = y)) +
geom_line() +
geom_area(
data = shaded_data_unif,
aes(x = x, y = y),
fill = "#ff8800",
alpha = 0.3
) +
geom_vline(xintercept = shade_up_to_x)
}
Reading paper:
Bai, J., So, K. C., Tang, C. S., Chen, X., & Wang, H. (2019). Coordinating supply and demand on an on-demand service platform with impatient customers. Manufacturing & Service Operations Management, 21(3), 556-570. https://doi.org/10.1287/msom.2018.0707
No focus on surge pricing for ease of modeling and cites customer dissatisfaction if price jumps in short time span.
Assumed \(p\) and \(w\) are known to customers and agents in advance.
Goal of service platform is to set \(p^*\) and \(w^*\) to maximize average profit.
Let \(\lambda\) be the customer request rate with max amount of service in a time period is \(\bar{\lambda}\)
Customer value per service unit is \(v\) (heterogeneous with CDF \(F(\cdot)\))
For a customer with valuation \(v\) and service request of \(D\) units the surplus is \((v - p)D\)
\(D\) is assumed to be independent of customer-type \(v\)
expected value of service units will be modeled as \(d = E(D)\)
Utility function of a customer can be modeled as \(U(v) = (v - p)d - cW_q\)
\(v\) is value per service unit
\(p\) is price per service unit
\(d\) is expected value of service units
\(c\) is cost of waiting
\(W_q\) is expected wait time for service
Realized customer request rate \(\lambda\) is the percent of \(\bar{\lambda}\) of customers who have \(U(v) \ge 0\) for a given \(p\)
aka \(\lambda = \text{Prob} \left\{ U(v) \ge 0 \right\} \cdot \bar{\lambda} \to \lambda = \text{Prob} \left\{ v \ge p + \frac{c}{d} W_q \right\} \cdot \bar{\lambda}\)
steps
\((v-p)d-cW_q \ge 0\)
\(vd-pd-cW_q \ge 0\)
\(vd \ge pd-cW_q\)
\(v \ge p - \frac{c}{d}W_q\)
This probability is labeled \(s\) (i.e. \(s = \text{Prob} \left\{ v \ge p + \frac{c}{d} W_q \right\}\) and \(\lambda = s \bar{\lambda}\)) and represents the desired service rate. We can manipulate \(p\) to decide what percentage of \(\bar{\lambda}\) to target
Since \(v \sim F(\cdot)\) then the price rate satisfies \(p = F^{-1} \left( 1 - \frac{\lambda}{\bar{\lambda}} \right) - \frac{c}{d} W_q\) (note this could be written \(p = F^{-1} \left( 1 - s \right) - \frac{c}{d} W_q\) to show service rate more directly)
Steps
\(\lambda = \text{Prob} \left\{ v \ge p + \frac{c}{d} W_q \right\} \cdot \bar{\lambda}\)
\(\lambda = \left[ 1 - F\left( p + \frac{c}{d} W_q \right) \right] \cdot \bar{\lambda}\)
\(\frac{\lambda}{\bar{\lambda}} = 1 - F\left( p + \frac{c}{d} W_q \right)\)
\(F\left( p + \frac{c}{d} W_q \right) = 1 - \frac{\lambda}{\bar{\lambda}}\)
\(p + \frac{c}{d} W_q = F^{-1}\left(1 - \frac{\lambda}{\bar{\lambda}}\right)\)
\(p = F^{-1}\left(1 - \frac{\lambda}{\bar{\lambda}}\right) - \frac{c}{d} W_q\) or \(p = F^{-1}\left(1 - s\right) - \frac{c}{d} W_q\)
Holding everything else constant
Price rate decreases as wait time ( \(W_q\) ) increases
Price rate decreases as unit wait cost ( \(c\) ) increases
Price rate increases as number of service units ( \(d\) ) increases
Price is prescribed as:
\[ p = F^{-1}\left(1 - s\right) - \frac{c}{d} W_q = F^{-1}\left(1 - \frac{\lambda}{\bar{\lambda}}\right) - \frac{c}{d} W_q \]
For given
# Price rate as a function of service level
price_rate <- function(lambda, lambda_bar, c, d, Wq, Finv = \(p) qunif(p = p, min = 0, max = 1)) {
s <- lambda / lambda_bar
Finv(1 - s) - (c / d) * Wq
}
Holding all other factors constant and just manipulating the input shown on x axis.
These relationships are the similar in behavior regardless of choice of value distribution family ( \(F(\cdot)\) ).
These notes don’t reflect how \(\lambda\) and \(W_q\) are related (i.e. more requests would lead to more wait time assuming the same supply).
ex_range <- 0:10
lambda_bar <- 10
lambda <- 2
c <- 2
d <- 2
Wq <- 2
# value dist
Finv <- \(p) qunif(p = p, min = 0, max = 1)
ex_lambda <- price_rate(ex_range, lambda_bar, ex_range, d, Wq, Finv)
ex_c <- price_rate(lambda, lambda_bar, ex_range, d, Wq, Finv)
ex_d <- price_rate(lambda, lambda_bar, c, ex_range, Wq, Finv)
ex_Wq <- price_rate(lambda, lambda_bar, c, d, ex_range, Finv)
func_behavior_df <- data.frame(
ex_range = ex_range,
lambda = ex_lambda,
c = ex_c,
d = ex_d,
Wq = ex_Wq
)
func_behavior_df |>
pivot_longer(-ex_range) |>
mutate(name = ifelse(name == "c", "Cost of waiting (c)", name)) |>
mutate(name = ifelse(name == "d", "Average service units (d)", name)) |>
mutate(name = ifelse(name == "Wq", "Wait time (Wq)", name)) |>
mutate(name = ifelse(name == "lambda", "Customer requests (lambda)", name)) |>
ggplot(aes(x = ex_range, y = value)) +
geom_line() +
# geom_hline(yintercept = 0, linetype = "dotted") +
facet_wrap(~name) +
scale_x_continuous(breaks = seq(0, 10, 2)) +
theme(
axis.text.y = element_blank(), axis.ticks.y = element_blank()
) +
labs(
y = "price", x = "",
title = bquote("Behavior of price with inputs; " ~
p == F^
{
-1
} * (1 - lambda / bar(lambda) - frac(c, d) * W[q])),
subtitle = bquote("Holding all other factors constant arbitrarily at 2 & " ~ bar(lambda) == 10)
)
The paper will use a uniform distribution \([0, 1]\) to simplify analysis but insights extend to other families
c <- 2
d <- 2
Wq <- 2
# value dist (unifo) params
Finv <- \(p) qunif(p = p, min = 0, max = 1)
lambda <- 2
lambda_bar <- 10
s <- lambda / lambda_bar
p1 <- plot_shaded_unif(a = 0, b = 1, p = s) +
labs(
title = sprintf(
"For service level %.0f%% -> price = %.2f ",
s * 100, price_rate(lambda, lambda_bar, c, d, Wq, Finv)
),
subtitle = sprintf("Given c = %.1f; d = %.1f, Wq = %.1f", c, d, Wq),
x = "v",
y = ""
) +
theme(
plot.title = element_text(size = 10),
plot.subtitle = element_text(size = 8),
axis.text.y = element_blank(),
axis.ticks.y = element_blank()
)
lambda <- 6
lambda_bar <- 10
s <- lambda / lambda_bar
p2 <- plot_shaded_unif(a = 0, b = 1, p = s) +
labs(
title = sprintf(
"For service level %.0f%% -> price = %.2f ",
s * 100, price_rate(lambda, lambda_bar, c, d, Wq, Finv)
),
subtitle = sprintf("Given c = %.1f; d = %.1f, Wq = %.1f", c, d, Wq),
x = "v",
y = ""
) +
theme(
plot.title = element_text(size = 10),
plot.subtitle = element_text(size = 8),
axis.text.y = element_blank(),
axis.ticks.y = element_blank()
)
p1 + p2
c <- 2
d <- 2
Wq <- 2
# value dist (normal) params
mu <- 5
sd <- 1
Finv <- \(p) qnorm(p = p, mean = mu, sd = sd)
lambda <- 3
lambda_bar <- 10
s <- lambda / lambda_bar
p1 <- plot_shaded_normal(mean = mu, sd = sd, p = s) +
labs(
title = sprintf(
"For service level %.0f%% -> price = %.2f ",
s * 100, price_rate(lambda, lambda_bar, c, d, Wq, Finv)
),
subtitle = sprintf("Given c = %.1f; d = %.1f, Wq = %.1f", c, d, Wq),
x = "v",
y = ""
) +
theme(
plot.title = element_text(size = 10),
plot.subtitle = element_text(size = 8),
axis.text.y = element_blank(),
axis.ticks.y = element_blank()
)
lambda <- 8
lambda_bar <- 10
s <- lambda / lambda_bar
p2 <- plot_shaded_normal(mean = mu, sd = sd, p = s) +
labs(
title = sprintf(
"For service level %.0f%% -> price = %.2f ",
s * 100, price_rate(lambda, lambda_bar, c, d, Wq, Finv)
),
subtitle = sprintf("Given c = %.1f; d = %.1f, Wq = %.1f", c, d, Wq),
x = "v",
y = ""
) +
theme(
plot.title = element_text(size = 10),
plot.subtitle = element_text(size = 8),
axis.text.y = element_blank(),
axis.ticks.y = element_blank()
)
p1 + p2
Wage is prescribed as:
\[ w = G^{-1}(\beta) \frac{k}{\lambda d} = G^{-1}(\frac{k}{K}) \frac{k}{\lambda d} \]
For given
# Wage rate as a function of provider participation rate
wage_rate <- function(k, K, lambda, d, Ginv = \(p) qunif(p = p, min = 0, max = 1)) {
Ginv(k / K) * (k / (lambda * d))
}
Holding all other factors constant and just manipulating the input shown on x axis.
These relationships are the similar in behavior regardless of choice of reservation rate distribution family ( \(G(\cdot)\) ).
ex_range <- 0:10
K <- 10
k <- 2
lambda <- 2
d <- 2
# reservation rate dist
Ginv <- \(p) qunif(p = p, min = 0, max = 1)
ex_k <- wage_rate(ex_range, K, lambda, d, Ginv)
ex_lambda <- wage_rate(k, K, ex_range, d, Ginv)
ex_d <- wage_rate(k, K, lambda, ex_range, Ginv)
func_behavior_df <- data.frame(
ex_range = ex_range,
k = ex_k,
lambda = ex_lambda,
d = ex_d
)
func_behavior_df |>
pivot_longer(-ex_range) |>
mutate(name = ifelse(name == "k", "Pariticpating providers (k)", name)) |>
mutate(name = ifelse(name == "d", "Average service units (d)", name)) |>
mutate(name = ifelse(name == "lambda", "Customer requests (lambda)", name)) |>
ggplot(aes(x = ex_range, y = value)) +
geom_line() +
# geom_hline(yintercept = 0, linetype = "dotted") +
facet_wrap(~name) +
scale_x_continuous(breaks = seq(0, 10, 2)) +
theme(
axis.text.y = element_blank(), axis.ticks.y = element_blank()
) +
labs(
y = "wage", x = "",
title = bquote("Behavior of wage with inputs; " ~
w == G^
{
-1
} * (frac(k, K)) * frac(k, lambda * d)),
subtitle = bquote("Holding all other factors constant arbitrarily at 2 & " ~ K == 10)
)
The paper will use a uniform distribution \([0, 1]\) to simplify analysis but insights extend to other families
lambda <- 2
d <- 2
# value dist (unifo) params
Ginv <- \(p) qunif(p = p, min = 0, max = 1)
k <- 2
K <- 10
beta <- k / K
p1 <- plot_shaded_unif(a = 0, b = 1, p = beta) +
labs(
title = sprintf(
"For partipication level %.0f%% -> wage = %.2f ",
beta * 100, wage_rate(k, K, lambda, d, Ginv)
),
subtitle = sprintf("Given lambda = %.1f; d = %.1f", lambda, d),
x = "v",
y = ""
) +
theme(
plot.title = element_text(size = 10),
plot.subtitle = element_text(size = 8),
axis.text.y = element_blank(),
axis.ticks.y = element_blank()
)
k <- 8
K <- 10
beta <- k / K
p2 <- plot_shaded_unif(a = 0, b = 1, p = beta) +
labs(
title = sprintf(
"For service level %.0f%% -> price = %.2f ",
beta * 100, wage_rate(k, K, lambda, d, Ginv)
),
subtitle = sprintf("Given lambda = %.1f; d = %.1f", lambda, d),
x = "v",
y = ""
) +
theme(
plot.title = element_text(size = 10),
plot.subtitle = element_text(size = 8),
axis.text.y = element_blank(),
axis.ticks.y = element_blank()
)
p1 + p2
Wage rate ( \(w\) ) will be a proportion ( \(\alpha\) ) of price rate ( \(p\) ) - \(w = \alpha p\)
Waiting time \(W_q\) will be based on \(M/M/k\) queue with arrival rate \(\lambda\) and service rate \(\frac{\mu}{d}\)
For simplicity in base model:
The distribution of \(v\) (so far referred to as \(F(\cdot)\) ) and \(r\) (so far referred to as \(G(\cdot)\) ) will both be modeled as uniform on range \([0, 1]\). This means we can drop out mention of them from current model.
\(p = F^{-1}\left(1 - \frac{\lambda}{\bar{\lambda}}\right) - \frac{c}{d} W_q\) becomes \(p = \left(1 - \frac{\lambda}{\bar{\lambda}}\right) - \frac{c}{d} W_q\)
\(w = G^{-1}(\frac{k}{K}) \frac{k}{\lambda d}\) becomes \(w = \frac{k}{K} \cdot \frac{k}{\lambda d} = \frac{k^2}{K \lambda d}\)
\(\max_{k, \lambda} \pi(k, \lambda) = \lambda d (p - w) = \lambda d \left[ \left( F^{-1}\left(1 - s\right) - \frac{c}{d} W_q \right) - \left( G^{-1}(\beta) \frac{k}{\lambda d} \right) \right]\) becomes \(\max_{k, \lambda} \pi(k, \lambda) = \lambda d (p - w) = \lambda d \left[ \left( \left(1 - s \right) - \frac{c}{d} W_q \right) - \left( \frac{k^2}{K \lambda d} \right) \right]\)
Wage rate is \(\alpha\) proportion of \(p\) : \(w = \alpha p \to \frac{w}{\alpha} = p\) ; avoid direct \(p\) formulation since it becomes way messier when incorporating Wait time formulation
\(\max_{k, \lambda} \pi(k, \lambda) = \lambda d (p - w) = \lambda d \left( \frac{w}{\alpha} - w \right)\)
Incorporating definition of \(w = \frac{k^2}{K \lambda d}\) we can write as \(\max_{k, \lambda} \pi(k, \lambda) = \lambda d \left( \frac{k^2}{K \lambda d \alpha} - \frac{k^2}{K \lambda d} \right) = \frac{k^2}{K \alpha} - \frac{k^2}{K} = \frac{k^2 - k^2 \alpha}{K \alpha} = \frac{k^2 (1 - \alpha)}{K \alpha}\)
Pretty clear to see we want to maximize number of participating providers \(k\)
This is subject to constraints
\(\frac{\lambda d}{k \mu} < 1\)
\(w = \alpha p\) which can be written as (still avoiding writing out \(W_q\) here; see def above)
\(\frac{k^2}{K \lambda d} = \alpha \left[ \left(1 - \frac{\lambda}{\bar{\lambda}}\right) - \frac{c}{d} W_q \right]\)
This implies \(k^2 = {K \lambda d}\alpha \left[ \left(1 - \frac{\lambda}{\bar{\lambda}}\right) - \frac{c}{d} W_q \right]\)
exact_wait_time_ <- function(k, lambda, mu, d) {
# $$
# W_q = \frac{1}{1 + \left( \frac{k! (1 - \rho)}{k^k \rho^k} \right) \sum_{i=0}^{k-1} \frac{k^i \rho^i}{i!}} \left[ \frac{\rho}{\lambda (1 - \rho)} \right]
# $$
# Calculate rho (system utilization)
rho <- (lambda * d) / (k * mu)
stopifnot(rho < 1)
# First fraction term
# ... Denominator term 1
## scalar to adjust probs of term 2
term1 <- (factorial(k) * (1 - rho)) / (k^k * rho^k)
# ... Denominator term 2
## each term in summation is prob of i servers busy
i <- 0:(k - 1)
term2 <- sum((k^i * rho^i) / factorial(i))
# Full formula for Wq
wait_time <- 1 / (1 + term1 * term2) * (rho / (lambda * (1 - rho)))
return(wait_time)
}
exact_wait_time <- Vectorize(exact_wait_time_, c("k", "lambda", "mu", "d"))
prob_queue_len <- function(i, k, lambda, mu, d) {
rho <- (lambda * d) / (k * mu)
n <- 0:(k - 1)
p_0 <- 1 / (sum(rho^n / factorial(n)) + (rho^k / factorial(k)) * (1 / (1 - rho)))
if (i == 0) {
p_i <- p_0
} else if (i < k) {
p_i <- (rho^i / factorial(i)) * p_0
} else {
p_i <- (rho^i / (factorial(k) * k^(i - k))) * p_0
}
return(p_i)
}
queue_len_probs <- function(k, lambda, mu, d) {
probs <- c()
i_vals <- 0:(k + 10)
for (i in i_vals) {
p_i <- prob_queue_len(i, k, lambda, mu, d)
probs[i + 1] <- p_i
}
return(data.frame(len = i_vals, prob = probs))
}
check_rho_constraint <- function(lambda, d, k, mu) {
rho <- (lambda * d) / (k * mu)
return(rho < 1)
}
check_fixed_payout_constraint <- function(
k, K, lambda, lambda_bar,
c, d,
mu,
alpha,
tolerance = 0.01) {
v1 <- k^2
Wq <- exact_wait_time(k, lambda, mu, d)
v2 <- K * lambda * d * alpha * ((1 - (lambda / lambda_bar)) - (c / d) * Wq)
return(isTRUE(all.equal(v1, v2, tolerance = tolerance)))
}
profit_function <- function(
k, K, lambda, lambda_bar,
c, d,
mu,
alpha,
Finv = \(p) qunif(p, min = 0, max = 1),
Ginv = \(p) qunif(p, min = 0, max = 1)) {
meets_rho_constraint <- check_rho_constraint(lambda, d, k, mu)
if (!meets_rho_constraint) {
return(list(
meets_rho_constraint = FALSE,
meets_fixed_payout_constraint = FALSE
))
}
meets_fixed_payout_constraint <- check_fixed_payout_constraint(k, K, lambda, lambda_bar, c, d, mu, alpha)
if (!meets_fixed_payout_constraint) {
return(list(
meets_rho_constraint = TRUE,
meets_fixed_payout_constraint = FALSE
))
}
Wq <- exact_wait_time(k, lambda, mu, d)
price <- price_rate(lambda, lambda_bar, c, d, Wq, Finv)
wage <- wage_rate(k, K, lambda, d, Ginv)
profit <- (k^2 * (1 - alpha)) / (K * alpha)
list(
k = k,
K = K,
lambda = lambda,
lambda_bar = lambda_bar,
c = c,
d = d,
mu = mu,
rho = (lambda * d) / (k * mu),
alpha = alpha,
Wq = Wq,
Lq = lambda * Wq,
queue_len_dist = queue_len_probs(k, lambda, mu, d),
price = price,
wage = wage,
profit = profit,
meets_rho_constraint = meets_rho_constraint,
meets_fixed_payout_constraint = meets_fixed_payout_constraint
)
}
display_results <- function(res) {
display_res <- res[c("k", "K", "lambda", "lambda_bar", "rho", "Wq", "price", "wage", "alpha", "profit")]
display_res |>
as.data.frame() |>
t() |>
as.data.frame() |>
rownames_to_column() |>
set_names(c("Param", "Value")) |>
mutate(Value = as.numeric(sprintf("%.3f", Value)))
}
Here is the results for row 1 of Table 1 to show my formulation matches.
res <- profit_function(k = 7, K = 50, lambda = 2.71, lambda_bar = 10, c = 1, d = 1, mu = 1, alpha = 0.5)
res$queue_len_dist |>
filter(prob > 0.01) |>
ggplot(aes(x = len, y = prob)) +
geom_bar(stat = "identity")
display_results(res)
Interestingly, and unreported by the paper. This is not a unique solution. Their Table 1 seems to prefer a solution with higher price and lower wait times.
res <- profit_function(k = 7, K = 50, lambda = 4.73, lambda_bar = 10, c = 1, d = 1, mu = 1, alpha = 0.5)
res$queue_len_dist |>
filter(prob > 0.01) |>
ggplot(aes(x = len, y = prob)) +
geom_bar(stat = "identity")
display_results(res)
My formulation matches paper results, but to save time I’ve used a lower precision. Because of this, my optimal values for lambda will be off from paper table by decimal amounts.
find_optimal_k_and_lambda <- function(
K, lambda_bar, c, d, mu, alpha,
Finv = \(p) qunif(p, min = 0, max = 1),
Ginv = \(p) qunif(p, min = 0, max = 1),
prefer_smaller_price = TRUE,
k_step = 1,
lambda_step = 0.01,
tied_profit_tolerance = 0.001) {
best_outcome <- NULL
possible_k <- seq(1, K, by = k_step)
# if (interactive()) {
# n_iter <- length(possible_k) * length(possible_lambda)
# pb <- txtProgressBar(min = 0, max = n_iter, style = 3)
# progress <- 0
# }
for (k in possible_k) {
possible_lambda <- seq(lambda_step, k, by = lambda_step)
for (lambda in possible_lambda) {
res <- profit_function(
k = k, K = K,
lambda = lambda, lambda_bar = lambda_bar,
c = c, d = d,
mu = mu,
alpha = alpha,
Finv = Finv, Ginv = Ginv
)
if (res$meets_rho_constraint & res$meets_fixed_payout_constraint) {
if (is.null(best_outcome)) {
best_outcome <- res
next
}
if (res$profit > best_outcome$profit) {
best_outcome <- res
next
}
# if profit is tied with best (within a tolerance)
if (isTRUE(all.equal(res$profit, best_outcome$profit, tolerance = tied_profit_tolerance))) {
if (prefer_smaller_price & res$price < best_outcome$price) {
best_outcome <- res
} else if (!prefer_smaller_price & res$price > best_outcome$price) {
best_outcome <- res
}
}
}
# if (interactive()) {
# progress <- progress + 1
# setTxtProgressBar(pb, progress)
# }
}
}
# if (interactive()) {
# close(pb)
# }
return(best_outcome)
}
Here is the results for row 1 of Table 1 to show my formulation can find the same solution, but with a lack of precision.
res <- find_optimal_k_and_lambda(
K = 50,
lambda_bar = 10,
c = 1,
d = 1,
mu = 1,
alpha = 0.5,
prefer_smaller_price = FALSE
)
display_results(res)
res <- find_optimal_k_and_lambda(
K = 50,
lambda_bar = 10,
c = 1,
d = 1,
mu = 1,
alpha = 0.5,
prefer_smaller_price = TRUE
)
display_results(res)
all_res <- list()
i <- 0
for (alpha in seq(0.1, 0.9, 0.1)) {
i <- i + 1
res <- find_optimal_k_and_lambda(
K = 50,
lambda_bar = 10,
c = 1,
d = 1,
mu = 1,
alpha = alpha,
prefer_smaller_price = FALSE
)
res$queue_len_dist <- NULL
all_res[[i]] <- res
}
res_df <- do.call(rbind, lapply(all_res, as.data.frame))
res_df
res_df |>
select(k, lambda, rho, alpha, Wq, price, wage, profit) |>
pivot_longer(-alpha) |>
ggplot(aes(x = alpha, y = value)) +
geom_line() +
geom_point() +
facet_wrap(~name, scales = "free_y")
all_res <- list()
i <- 0
for (alpha in seq(0.1, 0.9, 0.1)) {
i <- i + 1
res <- find_optimal_k_and_lambda(
K = 50,
lambda_bar = 10,
c = 1,
d = 1,
mu = 1,
alpha = alpha,
prefer_smaller_price = TRUE
)
res$queue_len_dist <- NULL
all_res[[i]] <- res
}
res_df <- do.call(rbind, lapply(all_res, as.data.frame))
res_df
res_df |>
select(k, lambda, rho, alpha, Wq, price, wage, profit) |>
pivot_longer(-alpha) |>
ggplot(aes(x = alpha, y = value)) +
geom_line() +
geom_point() +
facet_wrap(~name, scales = "free_y")
all_res <- list()
i <- 0
for (lambda_bar in seq(10, 100, 10)) {
i <- i + 1
res <- find_optimal_k_and_lambda(
K = 50,
lambda_bar = lambda_bar,
c = 1,
d = 1,
mu = 1,
alpha = 0.5,
prefer_smaller_price = FALSE
)
res$queue_len_dist <- NULL
all_res[[i]] <- res
}
res_df <- do.call(rbind, lapply(all_res, as.data.frame))
res_df
res_df |>
select(k, lambda, lambda_bar, rho, Wq, price, wage, profit) |>
pivot_longer(-lambda_bar) |>
ggplot(aes(x = lambda_bar, y = value)) +
geom_line() +
geom_point() +
facet_wrap(~name, scales = "free_y")
all_res <- list()
i <- 0
for (lambda_bar in seq(10, 100, 10)) {
i <- i + 1
res <- find_optimal_k_and_lambda(
K = 50,
lambda_bar = lambda_bar,
c = 1,
d = 1,
mu = 1,
alpha = 0.5,
prefer_smaller_price = TRUE
)
res$queue_len_dist <- NULL
all_res[[i]] <- res
}
res_df <- do.call(rbind, lapply(all_res, as.data.frame))
res_df
res_df |>
select(k, lambda, lambda_bar, rho, Wq, price, wage, profit) |>
pivot_longer(-lambda_bar) |>
ggplot(aes(x = lambda_bar, y = value)) +
geom_line() +
geom_point() +
facet_wrap(~name, scales = "free_y")
all_res <- list()
i <- 0
for (d in seq(1, 10, 2)) {
i <- i + 1
res <- find_optimal_k_and_lambda(
K = 50,
lambda_bar = 10,
c = 1,
d = d,
mu = 1,
alpha = 0.5,
prefer_smaller_price = FALSE
)
res$queue_len_dist <- NULL
all_res[[i]] <- res
}
res_df <- do.call(rbind, lapply(all_res, as.data.frame))
res_df
res_df |>
select(k, lambda, d, rho, Wq, price, wage, profit) |>
pivot_longer(-d) |>
ggplot(aes(x = d, y = value)) +
geom_line() +
geom_point() +
facet_wrap(~name, scales = "free_y")
all_res <- list()
i <- 0
for (d in seq(1, 10, 2)) {
i <- i + 1
res <- find_optimal_k_and_lambda(
K = 50,
lambda_bar = 10,
c = 1,
d = d,
mu = 1,
alpha = 0.5,
prefer_smaller_price = TRUE
)
res$queue_len_dist <- NULL
all_res[[i]] <- res
}
res_df <- do.call(rbind, lapply(all_res, as.data.frame))
res_df
res_df |>
select(k, lambda, d, rho, Wq, price, wage, profit) |>
pivot_longer(-d) |>
ggplot(aes(x = d, y = value)) +
geom_line() +
geom_point() +
facet_wrap(~name, scales = "free_y")
all_res <- list()
i <- 0
for (alpha in seq(0.1, 0.9, 0.1)) {
i <- i + 1
res <- find_optimal_k_and_lambda(
K = 50,
lambda_bar = 20,
c = 1,
d = 1,
mu = 1,
alpha = alpha,
k_step = 0.2,
lambda_step = 0.2,
prefer_smaller_price = FALSE
)
res$queue_len_dist <- NULL
all_res[[i]] <- res
}
res_df <- do.call(rbind, lapply(all_res, as.data.frame))
res_df
res_df |>
select(k, lambda, rho, alpha, Wq, price, wage, profit) |>
pivot_longer(-alpha) |>
ggplot(aes(x = alpha, y = value)) +
geom_line() +
geom_point() +
facet_wrap(~name, scales = "free_y")
all_res <- list()
i <- 0
for (lambda_bar in seq(10, 100, 10)) {
i <- i + 1
res <- find_optimal_k_and_lambda(
K = 50,
lambda_bar = lambda_bar,
c = 1,
d = 1,
mu = 1,
alpha = 0.5,
k_step = 0.2,
lambda_step = 0.2,
prefer_smaller_price = FALSE
)
res$queue_len_dist <- NULL
all_res[[i]] <- res
}
res_df <- do.call(rbind, lapply(all_res, as.data.frame))
res_df
res_df |>
select(k, lambda, lambda_bar, rho, Wq, price, wage, profit) |>
pivot_longer(-lambda_bar) |>
ggplot(aes(x = lambda_bar, y = value)) +
geom_line() +
geom_point() +
facet_wrap(~name, scales = "free_y")
all_res <- list()
i <- 0
for (K in seq(10, 100, 10)) {
i <- i + 1
res <- find_optimal_k_and_lambda(
K = K,
lambda_bar = 20,
c = 1,
d = 1,
mu = 1,
alpha = 0.5,
k_step = 0.2,
lambda_step = 0.2,
prefer_smaller_price = FALSE
)
res$queue_len_dist <- NULL
all_res[[i]] <- res
}
res_df <- do.call(rbind, lapply(all_res, as.data.frame))
res_df
res_df |>
select(k, lambda, K, rho, Wq, price, wage, profit) |>
pivot_longer(-K) |>
ggplot(aes(x = K, y = value)) +
geom_line() +
geom_point() +
facet_wrap(~name, scales = "free_y")
all_res <- list()
i <- 0
for (K in seq(10, 100, 20)) {
for (lambda_bar in seq(10, 100, 20)) {
i <- i + 1
res <- find_optimal_k_and_lambda(
K = K,
lambda_bar = lambda_bar,
c = 1,
d = 1,
mu = 1,
alpha = 0.5,
k_step = 1,
lambda_step = 1,
prefer_smaller_price = FALSE
)
res$queue_len_dist <- NULL
all_res[[i]] <- res
}
}
res_df <- do.call(rbind, lapply(all_res, as.data.frame))
res_df
plot_cols <- c("k", "lambda", "rho", "Wq", "price", "wage", "profit")
res_df |>
mutate_at(all_of(plot_cols), scale) |>
select(K, lambda_bar, k, lambda, rho, Wq, price, wage, profit) |>
pivot_longer(all_of(plot_cols)) |>
rename(`z-score` = value) |>
ggplot(aes(x = K, y = lambda_bar, fill = `z-score`)) +
geom_tile() +
facet_wrap(~name)
## Warning: Using `all_of()` outside of a selecting function was deprecated in tidyselect
## 1.2.0.
## ℹ See details at
## <https://tidyselect.r-lib.org/reference/faq-selection-context.html>
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
# ... to be continued
To approximate and give tractable analytical results can instead be approximated with \(W_q = \frac{ \rho ^ \sqrt{2 (k + 1)} }{ \lambda (1 - \rho) }\)
See below for how this approximation compares to more exact, but more unwieldy formulation.
The paper actually goes further from exact calc and rewrites to be \(W_q = \frac{ \rho ^ \sqrt{2 (n + 1)} }{ \lambda (1 - \rho) }\). And \(n\) is set to be \(k^*\). This means for \(k^*\), \(W_q\) is well approximated, but for values of \(k < k^*\) the \(W_q\) will be over estimated and for values of \(k > k^*\), \(W_q\) will be underestimated.
approx_wait_time_ <- function(k, lambda, mu, d) {
# $$
# W_q = \frac{ \rho ^ \sqrt{2 (k + 1)} }{ \lambda (1 - \rho) }
# $$
# Calculate rho (system utilization)
rho <- (lambda * d) / (k * mu)
# numerator shows traffic intensity - for large k, less sensitive to traffic
# denominator combines arrival rate and idle capacity
# to show traffic that can be handled without queueing
wait_time <- (rho^sqrt(2 * (k + 1))) / (lambda * (1 - rho))
return(wait_time)
}
approx_wait_time <- Vectorize(approx_wait_time_, c("k", "lambda", "mu", "d"))
k <- 1:50
lambda <- 4
mu <- 16
d <- 2
wait_time_compare <- data.frame(
k = k,
Exact = exact_wait_time(k, lambda, mu, d),
Approx = approx_wait_time(k, lambda, mu, d)
) |>
pivot_longer(-k, names_to = "Formula") |>
mutate(`log(Wait time)` = log(value)) |>
rename(`Wait time` = value) |>
pivot_longer(contains("Wait"), names_to = "trans") |>
mutate(trans = forcats::fct_rev(trans), Formula = forcats::fct_rev(Formula))
ggplot(wait_time_compare, aes(x = k, y = value, color = Formula, linetype = Formula)) +
geom_line() +
facet_wrap(~trans, scales = "free_y") +
labs(
title = "Comparison of Approx and Exact M/M/k wait time formulas",
subtitle = bquote(lambda == .(lambda) ~ ", " ~
mu == .(mu) ~ ", " ~
d == .(d) ~ " --> " ~
rho == frac(lambda * d, k * mu) ~ ">=" ~ 1)
)